##==============================================================================
## This script replicates the main results presented in "Anchoring vignettes as
## a diagnostic tool for cross-national (in)comparability of survey measures: 
## The case of voters left-right self-placement."             
##                                                                               
## Nick Lin and Seonghui Lee
## British Journal of Political Science, accepted April 2023
##==============================================================================

rm(list=ls())
library(foreign) 
library(stringr)    ## ver 1.4.0
library(reshape2)   ## ver 0.8.8
library(dplyr)      ## ver 1.0.10
library(doBy)       ## ver 4.6.11
library(ggplot2)    ## ver 3.3.5
library(lme4)       ## ver 1.1-27.1 
options(scipen = 999)

################################################################################
##  Table 1 
################################################################################

###======================================
### Political Efficacy (King et al. 2004)
###======================================

load("Replication_King2004.Rda")

#### Restructure the data so that the unit is individual-vignette. 
ind.long.eff <- melt(king.data, id.vars = c("respID", "country"))

## Perform the Model
mod.rd <- lmer(value ~ (1|country) + (1|variable) + (1|variable:country), data = ind.long.eff, REML = TRUE)
summary(mod.rd)
re_dat = as.data.frame(VarCorr(mod.rd))
var.explained <- (re_dat[1,4] + re_dat[3,4])  / (re_dat[1,4] + re_dat[2,4] + re_dat[3,4])
var.explained #0.903


###======================================
### Democracy (Bratton 2010)
###======================================

load("Replication_Bratton2010.Rda")

## Restructure the data so that the unit is individual-vignette. 
ind.long.dem <- melt(dem.data, id.vars = c("respno", "country"))

## Perform the Model
mod.rd <- lmer(value ~ (1|country) + (1|variable) + (1|variable:country), data = ind.long.dem, REML = TRUE)
summary(mod.rd)
re_dat = as.data.frame(VarCorr(mod.rd))
var.explained <- (re_dat[1,4] + re_dat[2,4])  / (re_dat[1,4] + re_dat[2,4] + re_dat[3,4])
var.explained #0.054


###======================================
### Expert's LR Placement (Bakker et al. 2014)
###======================================

load("Replication_Bakker2014.Rda")

#### Restructure the data so that the unit is individual-vignette. 
ind.long.expert <- melt(exp.dat, id.vars = c("respID", "country"))

## Perform the Model
mod.rd <- lmer(value ~ (1|country) + (1|variable) + (1|variable:country), data = ind.long.expert, REML = TRUE)
summary(mod.rd)
re_dat = as.data.frame(VarCorr(mod.rd))
var.explained <- (re_dat[1,4] + re_dat[2,4])  / (re_dat[1,4] + re_dat[2,4] + re_dat[3,4])
var.explained #0.011


###======================================
### Political Interest (Lee et al. 2015)
###======================================

load("Replication_LLS2015.Rda")

#### Restructure the data so that the unit is individual-vignette. 
ind.long.int.15 <- melt(int15, id.vars = c("resID", "country"))

## Perform the Model
mod.rd <- lmer(value ~ (1|country) + (1|variable) + (1|variable:country), data = ind.long.int.15, REML = TRUE)
summary(mod.rd)
re_dat = as.data.frame(VarCorr(mod.rd))
var.explained <- (re_dat[1,4] + re_dat[3,4])  / (re_dat[1,4] + re_dat[2,4] + re_dat[3,4])
var.explained #0.069


###======================================
### Political Interest (Lee et al. 2016)
###======================================

load("Replication_LLS2016.Rda")

#### Restructure the data so that the unit is individual-vignette. 
ind.long.int.16 <- melt(int16, id.vars = c("respid", "country"))

## Perform the Model
mod.rd <- lmer(value ~ (1|country) + (1|variable) + (1|variable:country), data = ind.long.int.16, REML = TRUE)
summary(mod.rd)
re_dat = as.data.frame(VarCorr(mod.rd))
var.explained <- (re_dat[1,4] + re_dat[2,4])  / (re_dat[1,4] + re_dat[2,4] + re_dat[3,4])
var.explained #0.033


###======================================
### LR Placement
###======================================

load("Replication_LR_Vign.Rda")

#### Restructure the data so that the unit is individual-vignette. 
long.lr.20 <- melt(lr.data.20, id.vars = c("respID", "country"))

## Perform the Randon-Effect Model 
mod.rd <- lmer(value ~ (1|country) + (1|variable) + (1|variable:country), data = long.lr.20, REML = TRUE)
summary(mod.rd)
re_dat = as.data.frame(VarCorr(mod.rd))
var.explained <- (re_dat[1,4] + re_dat[2,4])  / (re_dat[1,4] + re_dat[2,4] + re_dat[3,4])
var.explained #0.112



################################################################################
##  Figure 2 
##  Note: the codes below use the data we have generated through bootstrapping 
##  process as the process demands lengthy computing time. The codes for detailed
##  bootstrapping process are found next section of the script.
################################################################################

## For codes producing the following CNDIF data sets, please see lines 174-620.
  load("CNDIF_LR_Vign.Rda")
  load("CNDIF_LLS16.Rda")
  load("CNDIF_Bratton10.Rda")
  load("CNDIF_Bakker14.Rda")
##
  
levels(exp.dif.CI$country)[levels(exp.dif.CI$country)=="United Kingdom"] <- "UK" 
exp.dif.CI$study <- "Bakker et al. (2014)" 
lr.dif.CI.20$study <- "Left-Right Placement (Our Study)" 
int.dif.CI$study <- "Lee et al. (2016)" 
dem.dif.CI$study <- "Bratton (2010)" 

all.data <- rbind.data.frame(exp.dif.CI, lr.dif.CI.20, int.dif.CI, dem.dif.CI)
all.data$country <- factor(all.data$country, 
                           levels=c("Argentina", "Austria", "Belgium", "Bulgaria", "Canada", "Chile", "China", "Czech Rep.", "Denmark", "Estonia", "Finland",
                                    "France", "Germany", "Greece", "Hungary", "Ireland", "Italy", "Japan", "Korea", "Latvia", "Lithuania", "Mexico",
                                    "Netherlands", "New Zealand", "Norway", "Poland", "Portugal", "Romania", "Slovakia", "Slovenia", 
                                    "Spain", "Sweden", "Switzerland", "UK", "US", "Benin", "Botswana", "Burkina Faso", "Cape Verde", "Ghana", "Kenya", "Lesotho", "Liberia", "Madagascar",
                                    "Malawi","Mali", "Mozambique", "Namibia", "Nigeria", "Senegal", "South Africa", "Tanzania", "Uganda",
                                    "Zambia", "Zimbabwe"))
all.data$country <- factor(all.data$country, levels=rev(levels(all.data$country)))

plot.dif <- ggplot(all.data, aes(x = mean, xmin = lwb.nr95, xmax = upb.nr95, y = factor(country))) +
  geom_point() + theme_bw() + 
  geom_segment(aes(x = lwb.nr95, xend = upb.nr95, y = factor(country), yend=factor(country))) +
  facet_wrap(~study, ncol=2, scales = "free") +
  labs(x = "Country Specific CN-DIF Score", y = "") + 
  theme(axis.title.x = element_text(size = 14, vjust = .3))+ 
  theme(axis.text=element_text(size=12)) +
  scale_x_continuous(limits =c(0, .5), breaks=c(0, .25, .5)) +  
  theme(strip.text.x = element_text(size = 14, colour = "black", face="bold"))

ggsave(file="Fig2.tiff", plot.dif, height=9, width=12, units="in", dpi=600)



################################################################################
##  Figure 2 
##  Note: this section includes codes for detailed bootstrapping process that 
##  generates data for Figure 2. For a simpler replication of Figure 2, see the 
##  codes above.
################################################################################

##############################################
### Upper-Left Panel: Bakker et al. (2014) ###
##############################################

load("Replication_Bakker2014.Rda")

## Generate Every Possible Combination of Countries

ctrylist <- c("Austria", "Belgium", "Bulgaria", "Czech Rep.", "Denmark", "Estonia", "Finland",
              "France", "Germany", "Greece", "Hungary", "Ireland", "Italy", "Latvia", "Lithuania",
              "Netherlands", "Norway", "Poland", "Portugal", "Romania", "Slovakia", "Slovenia", 
              "Spain", "Sweden", "Switzerland", "United Kingdom")

combnList <- list()
for (i in 1:length(ctrylist)) {
  combnList[[i]] <- combn(x = ctrylist, m = i)
}

##### Prepare the boostrap sample
boot.country <- list()
boot.sample  <- list()
temp.dat <- list()
agg.dif.dat <- list()
ctry.dif.dat <- list()
## Bootstrap begins here
set.seed(77005)
for (b in 1:1000) { # b = number of boostrapped sample
  
  for (c in 1:length(ctrylist)) { ## This loop creates one boostrapped sample of all countries
    ctry <- subset(exp.dat, country==ctrylist[c]) ## sample one country at a time
    boot.country[[c]] <- sample_frac(ctry, size = 1, replace = TRUE) 
  }
  ## Collapse the sampled data to the country level, and store it. 
  ctry.sample <- do.call("rbind", boot.country) 
  boot.sample[[b]] <- melt(ctry.sample, id.vars = c("respID", "country"))
  
  ## Now perform a set of multilevel models for each country pair.
  ## We focus on country pair as this includes the most similar and most different combination.
  k=2 # country number = 2 in a combination
  
  for (i in 1:ncol(combnList[[k]])) { # i = the ith combination in comnList k
    ctryData <- list()
    dif <- NULL
    for (j in 1:k) {  # j = jth row in each combination
      ctryData[[j]]<- subset(boot.sample[[b]], country==combnList[[k]][j,i])
    }
    dat <- do.call("rbind", ctryData)
    mod <- lmer(value ~ (1|country) + (1|variable) + (1|variable:country), data = dat, REML = TRUE)
    summary(mod)
    # Generate the cross-national DIF 
    re_dat = as.data.frame(VarCorr(mod))
    dif <- (re_dat[1,4] + re_dat[3,4])  / (re_dat[1,4] + re_dat[2,4] + re_dat[3,4])
    temp.dat[[i]] <- cbind(dif, paste(combnList[[k]][,i], collapse=", "))
  }
  
  dif.dat <- do.call("rbind", temp.dat)
  colnames(dif.dat) <- c("DIF", "Countries")
  agg.dif.dat[[b]] <- data.frame(dif.dat)
  agg.dif.dat[[b]]$DIF.num <- as.numeric(agg.dif.dat[[b]]$DIF)
  
  ## The Country Specific DIF Index 
  dat <- list()
  for (c in 1:length(ctrylist)) {
    ctry = ctrylist[c]
    subset<- grep(ctry, agg.dif.dat[[b]]$Countries) 
    
    temp.dat <- list()
    
    for (q in 1:length(subset)) {  # j = jth row in each combination
      h = subset[q]
      temp.dat[[q]]<- agg.dif.dat[[b]][h,]
    }
    dat[[c]] <- do.call("rbind", temp.dat)
    dat[[c]]$country <- ctry
  }
  exp.agg.dif.ctry <- do.call("rbind", dat)
  
  ctry.dif.dat[[b]] <- summaryBy(DIF.num ~ country, data=exp.agg.dif.ctry, FUN=mean, na.rm=TRUE)
}

#### Boostrapped CN-DIF and CIs for each country    
boot.ctry.dif <- do.call("rbind", ctry.dif.dat)

boot.mdn <- NA
boot.upb.nr95 <- NA
boot.lwb.nr95 <- NA
boot.upb.pc95 <- NA
boot.lwb.pc95 <- NA
boot.mean<-NA
boot.country <-NA
for (c in 1:length(ctrylist)) { 
  ctry <- subset(boot.ctry.dif, country==ctrylist[c])
  boot.mean[c] <- mean(ctry$DIF.num.mean)
  se <- sd(ctry$DIF.num.mean)
  boot.lwb.nr95[c] <- mean(ctry$DIF.num.mean) - 1.96*se
  boot.upb.nr95[c] <- mean(ctry$DIF.num.mean) + 1.96*se
  
  boot.mdn[c]  <- quantile(ctry$DIF.num.mean, probs = .5)
  boot.lwb.pc95[c]  <- quantile(ctry$DIF.num.mean, probs = .025)
  boot.upb.pc95[c]  <- quantile(ctry$DIF.num.mean, probs = .975)
  boot.country[c] <- ctrylist[c]
}

exp.dif.CI <- cbind.data.frame(boot.lwb.nr95, boot.lwb.pc95, boot.mdn, boot.mean, boot.upb.nr95, boot.upb.pc95, boot.country)
colnames(exp.dif.CI) <- c("lwb.nr95", "lwb.pc95", "mdn", "mean", "upb.nr95", "upb.pc95", "country")
exp.dif.CI$study <- rep("Expert LR Placement", length(exp.dif.CI$country))
exp.dif.CI$lwb.nr95[exp.dif.CI$lwb.nr95 < 0 ] <- 0

## Save the bootstrapped data
save(exp.dif.CI, file="CNDIF_Bakker14.Rda")


###======================================
### Upper-Right Panel: Bratton (2010) 
###======================================

load("Replication_Bratton2010.Rda")

## Generate Every Possible Combination of Countries

ctrylist <- c("Benin", "Botswana", "Burkina Faso", "Cape Verde", "Ghana", "Kenya", "Lesotho", "Liberia", "Madagascar",
              "Malawi","Mali", "Mozambique", "Namibia", "Nigeria", "Senegal", "South Africa", "Tanzania", "Uganda",
              "Zambia", "Zimbabwe")

combnList <- list()
for (i in 1:length(ctrylist)) {
  combnList[[i]] <- combn(x = ctrylist, m = i)
}

##### Prepare the boostrap sample
boot.country <- list()
boot.sample  <- list()
temp.dat <- list()
agg.dif.dat <- list()
ctry.dif.dat <- list()
## Bootstrap begins here
set.seed(77005)
for (b in 1:1000) { # b = number of boostrapped sample

  for (c in 1:length(ctrylist)) { ## This loop creates one boostrapped sample of all countries
    ctry <- subset(dem.data, country==ctrylist[c]) ## sample one country at a time
    boot.country[[c]] <- sample_frac(ctry, size = 1, replace = TRUE) 
  }
  ## Collapse the sampled data to the country level, and store it. 
  ctry.sample <- do.call("rbind", boot.country) 
  boot.sample[[b]] <- melt(ctry.sample, id.vars = c("respno", "country"))
  
  ## Now perform a set of multilevel models for each country pair.
  ## We focus on country pair as this includes the most similar and most different combination.
  k=2 # country number = 2 in a combination
  
  for (i in 1:ncol(combnList[[k]])) { # i = the ith combination in comnList k
    ctryData <- list()
    dif <- NULL
    for (j in 1:k) {  # j = jth row in each combination
      ctryData[[j]]<- subset(boot.sample[[b]], country==combnList[[k]][j,i])
    }
    dat <- do.call("rbind", ctryData)
    mod <- lmer(value ~ (1|country) + (1|variable) + (1|variable:country), data = dat, REML = TRUE)
    summary(mod)
    # Generate the cross-national DIF 
    re_dat = as.data.frame(VarCorr(mod))
    dif <- (re_dat[1,4] + re_dat[3,4])  / (re_dat[1,4] + re_dat[2,4] + re_dat[3,4])
    temp.dat[[i]] <- cbind(dif, paste(combnList[[k]][,i], collapse=", "))
  }
  
  dif.dat <- do.call("rbind", temp.dat)
  colnames(dif.dat) <- c("DIF", "Countries")
  agg.dif.dat[[b]] <- data.frame(dif.dat)
  agg.dif.dat[[b]]$DIF.num <- as.numeric(agg.dif.dat[[b]]$DIF)
  
  ## The Country Specific DIF Index 
  dat <- list()
  for (c in 1:length(ctrylist)) {
    ctry = ctrylist[c]
    subset<- grep(ctry, agg.dif.dat[[b]]$Countries) 
    
    temp.dat <- list()
    
    for (q in 1:length(subset)) {  # j = jth row in each combination
      h = subset[q]
      temp.dat[[q]]<- agg.dif.dat[[b]][h,]
    }
    dat[[c]] <- do.call("rbind", temp.dat)
    dat[[c]]$country <- ctry
  }
  dem.agg.dif.ctry <- do.call("rbind", dat)
  
  ctry.dif.dat[[b]] <- summaryBy(DIF.num ~ country, data=dem.agg.dif.ctry, FUN=mean, na.rm=TRUE)
}

#### Boostrapped CN-DIF and CIs for each country    
boot.ctry.dif <- do.call("rbind", ctry.dif.dat)

boot.mdn <- NA
boot.upb.nr95 <- NA
boot.lwb.nr95 <- NA
boot.upb.pc95 <- NA
boot.lwb.pc95 <- NA
boot.mean<-NA
boot.country <-NA
for (c in 1:length(ctrylist)) { 
  ctry <- subset(boot.ctry.dif, country==ctrylist[c])
  boot.mean[c] <- mean(ctry$DIF.num.mean)
  se <- sd(ctry$DIF.num.mean)
  boot.lwb.nr95[c] <- mean(ctry$DIF.num.mean) - 1.96*se
  boot.upb.nr95[c] <- mean(ctry$DIF.num.mean) + 1.96*se
  
  boot.mdn[c]  <- quantile(ctry$DIF.num.mean, probs = .5)
  boot.lwb.pc95[c]  <- quantile(ctry$DIF.num.mean, probs = .025)
  boot.upb.pc95[c]  <- quantile(ctry$DIF.num.mean, probs = .975)
  boot.country[c] <- ctrylist[c]
}

dem.dif.CI <- cbind.data.frame(boot.lwb.nr95, boot.lwb.pc95, boot.mdn, boot.mean, boot.upb.nr95, boot.upb.pc95, boot.country)
colnames(dem.dif.CI) <- c("lwb.nr95", "lwb.pc95", "mdn", "mean", "upb.nr95", "upb.pc95", "country")
dem.dif.CI$study <- rep("Democratic Assessment", length(dem.dif.CI$country))

## Save the bootstrapped data
save(dem.dif.CI, file="CNDIF_Bratton10.Rda")


###======================================
### Bottom-Left Panel: Lee et al. (2016) 
###======================================

load("Replication_LLS2016.Rda")

## Generate Every Possible Combination of Countries

ctrylist <- c("Argentina", "Canada", "Chile", "China", "Denmark", "Germany", "Japan",
              "Korea", "Mexico", "New Zealand", "Poland", "US")

combnList <- list()
for (i in 1:length(ctrylist)) {
  combnList[[i]] <- combn(x = ctrylist, m = i)
}

##### Prepare the boostrap sample
boot.country <- list()
boot.sample  <- list()
temp.dat <- list()
agg.dif.dat <- list()
ctry.dif.dat <- list()
## Bootstrap begins here
set.seed(77005)
for (b in 1:1000) { # b = number of boostrapped sample
  
  for (c in 1:length(ctrylist)) { ## This loop creates one boostrapped sample of all countries
    ctry <- subset(int16, country==ctrylist[c]) ## sample one country at a time
    boot.country[[c]] <- sample_frac(ctry, size = 1, replace = TRUE) 
  }
  ## Collapse the sampled data to the country level, and store it. 
  ctry.sample <- do.call("rbind", boot.country) 
  boot.sample[[b]] <- melt(ctry.sample, id.vars = c("respid", "country"))
  
  ## Now perform a set of multilevel models for each country pair.
  ## We focus on country pair as this includes the most similar and most different combination.
  k=2 # country number = 2 in a combination
  
  for (i in 1:ncol(combnList[[k]])) { # i = the ith combination in comnList k
    ctryData <- list()
    dif <- NULL
    for (j in 1:k) {  # j = jth row in each combination
      ctryData[[j]]<- subset(boot.sample[[b]], country==combnList[[k]][j,i])
    }
    dat <- do.call("rbind", ctryData)
    mod <- lmer(value ~ (1|country) + (1|variable) + (1|variable:country), data = dat, REML = TRUE)
    summary(mod)
    # Generate the cross-national DIF 
    re_dat = as.data.frame(VarCorr(mod))
    dif <- (re_dat[1,4] + re_dat[3,4])  / (re_dat[1,4] + re_dat[2,4] + re_dat[3,4])
    temp.dat[[i]] <- cbind(dif, paste(combnList[[k]][,i], collapse=", "))
  }
  
  dif.dat <- do.call("rbind", temp.dat)
  colnames(dif.dat) <- c("DIF", "Countries")
  agg.dif.dat[[b]] <- data.frame(dif.dat)
  #agg.dif.dat[[b]]$DIF.num <- as.numeric(levels(agg.dif.dat[[b]]$DIF))[agg.dif.dat[[b]]$DIF]
  agg.dif.dat[[b]]$DIF.num <- as.numeric(agg.dif.dat[[b]]$DIF)
  
  ## The Country Specific DIF Index 
  dat <- list()
  for (c in 1:length(ctrylist)) {
    ctry = ctrylist[c]
    subset<- grep(ctry, agg.dif.dat[[b]]$Countries) 
    
    temp.dat <- list()
    
    for (q in 1:length(subset)) {  # j = jth row in each combination
      h = subset[q]
      temp.dat[[q]]<- agg.dif.dat[[b]][h,]
    }
    dat[[c]] <- do.call("rbind", temp.dat)
    dat[[c]]$country <- ctry
  }
  int.agg.dif.ctry <- do.call("rbind", dat)
  
  ctry.dif.dat[[b]] <- summaryBy(DIF.num ~ country, data=int.agg.dif.ctry, FUN=mean, na.rm=TRUE)
}

#### Boostrapped CN-DIF and CIs for each country    
boot.ctry.dif <- do.call("rbind", ctry.dif.dat)

boot.mdn <- NA
boot.upb.nr95 <- NA
boot.lwb.nr95 <- NA
boot.upb.pc95 <- NA
boot.lwb.pc95 <- NA
boot.mean<-NA
boot.country <-NA
for (c in 1:length(ctrylist)) { 
  ctry <- subset(boot.ctry.dif, country==ctrylist[c])
  boot.mean[c] <- mean(ctry$DIF.num.mean)
  se <- sd(ctry$DIF.num.mean)
  boot.lwb.nr95[c] <- mean(ctry$DIF.num.mean) - 1.96*se
  boot.upb.nr95[c] <- mean(ctry$DIF.num.mean) + 1.96*se
  
  boot.mdn[c]  <- quantile(ctry$DIF.num.mean, probs = .5)
  boot.lwb.pc95[c]  <- quantile(ctry$DIF.num.mean, probs = .025)
  boot.upb.pc95[c]  <- quantile(ctry$DIF.num.mean, probs = .975)
  boot.country[c] <- ctrylist[c]
}

int.dif.CI <- cbind.data.frame(boot.lwb.nr95, boot.lwb.pc95, boot.mdn, boot.mean, boot.upb.nr95, boot.upb.pc95, boot.country)
colnames(int.dif.CI) <- c("lwb.nr95", "lwb.pc95", "mdn", "mean", "upb.nr95", "upb.pc95", "country")
int.dif.CI$study <- rep("Political Interest", length(int.dif.CI$country))

## Save the bootstrapped data
save(int.dif.CI, file="CNDIF_LLS16.Rda")


###======================================
### Bottom-Left Panel: Left-Right Placement (Our Study) 
###======================================

load("Replication_LR_Vign.Rda")

## Generate Every Possible Combination of Countries

ctrylist <- c("France", "Germany", "Hungary", "Italy", "Netherlands", "Poland", "Spain", "Sweden", "UK")

combnList <- list()
for (i in 1:length(ctrylist)) {
  combnList[[i]] <- combn(x = ctrylist, m = i)
}

##### Prepare the boostrap sample
boot.country <- list()
boot.sample  <- list()
temp.dat <- list()
agg.dif.dat <- list()
ctry.dif.dat <- list()

## Bootstrap begins here
set.seed(77005)
for (b in 1:1000) { # b = number of boostrapped sample
  
  for (c in 1:length(ctrylist)) { ## This loop creates one boostrapped sample of all countries
    ctry <- subset(lr.data.20, country==ctrylist[c]) ## sample one country at a time
    boot.country[[c]] <- sample_frac(ctry, size = 1, replace = TRUE) 
  }
  ## Collapse the sampled data to the country level, and store it. 
  ctry.sample <- do.call("rbind", boot.country) 
  boot.sample[[b]] <- melt(ctry.sample, id.vars = c("respID", "country"))
  
  ## Now perform a set of multilevel models for each country pair.
  ## We focus on country pair as this includes the most similar and most different combination.
  k=2 # country number = 2 in a combination
  
  for (i in 1:ncol(combnList[[k]])) { # i = the ith combination in comnList k
    ctryData <- list()
    dif <- NULL
    for (j in 1:k) {  # j = jth row in each combination
      ctryData[[j]]<- subset(boot.sample[[b]], country==combnList[[k]][j,i])
    }
    dat <- do.call("rbind", ctryData)
    mod <- lmer(value ~ (1|country) + (1|variable) + (1|variable:country), data = dat, REML = TRUE)
    summary(mod)
    # Generate the cross-national DIF 
    re_dat = as.data.frame(VarCorr(mod))
    dif <- (re_dat[1,4] + re_dat[3,4])  / (re_dat[1,4] + re_dat[2,4] + re_dat[3,4])
    temp.dat[[i]] <- cbind(dif, paste(combnList[[k]][,i], collapse=", "))
  }
  
  dif.dat <- do.call("rbind", temp.dat)
  colnames(dif.dat) <- c("DIF", "Countries")
  agg.dif.dat[[b]] <- data.frame(dif.dat)
  #agg.dif.dat[[b]]$DIF.num <- as.numeric(levels(agg.dif.dat[[b]]$DIF))[agg.dif.dat[[b]]$DIF]
  agg.dif.dat[[b]]$DIF.num <- as.numeric(agg.dif.dat[[b]]$DIF)
  
  ## The Country Specific DIF Index 
  dat <- list()
  for (c in 1:length(ctrylist)) {
    ctry = ctrylist[c]
    subset<- grep(ctry, agg.dif.dat[[b]]$Countries) 
    
    temp.dat <- list()
    
    for (q in 1:length(subset)) {  # j = jth row in each combination
      h = subset[q]
      temp.dat[[q]]<- agg.dif.dat[[b]][h,]
    }
    dat[[c]] <- do.call("rbind", temp.dat)
    dat[[c]]$country <- ctry
  }
  lr.agg.dif.ctry <- do.call("rbind", dat)
  
  ctry.dif.dat[[b]] <- summaryBy(DIF.num ~ country, data=lr.agg.dif.ctry, FUN=mean, na.rm=TRUE)
}

#### Boostrapped CN-DIF and CIs for each country    
boot.ctry.dif <- do.call("rbind", ctry.dif.dat)

boot.mdn <- NA
boot.upb.nr95 <- NA
boot.lwb.nr95 <- NA
boot.upb.pc95 <- NA
boot.lwb.pc95 <- NA
boot.mean<-NA
boot.country <-NA

for (c in 1:length(ctrylist)) { 
  ctry <- subset(boot.ctry.dif, country==ctrylist[c])
  boot.mean[c] <- mean(ctry$DIF.num.mean)
  se <- sd(ctry$DIF.num.mean)
  boot.lwb.nr95[c] <- mean(ctry$DIF.num.mean) - 1.96*se
  boot.upb.nr95[c] <- mean(ctry$DIF.num.mean) + 1.96*se
  
  boot.mdn[c]  <- quantile(ctry$DIF.num.mean, probs = .5)
  boot.lwb.pc95[c]  <- quantile(ctry$DIF.num.mean, probs = .025)
  boot.upb.pc95[c]  <- quantile(ctry$DIF.num.mean, probs = .975)
  boot.country[c] <- ctrylist[c]
}

lr.dif.CI.20 <- cbind.data.frame(boot.lwb.nr95, boot.lwb.pc95, boot.mdn, boot.mean, boot.upb.nr95, boot.upb.pc95, boot.country)
colnames(lr.dif.CI.20) <- c("lwb.nr95", "lwb.pc95", "mdn", "mean", "upb.nr95", "upb.pc95", "country")
lr.dif.CI.20$study <- rep("LR Survey (Vigns go first)", length(lr.dif.CI.20$country))

## Save the bootstrapped data
save(lr.dif.CI.20, file="CNDIF_LR_Vign.Rda")

############################################################################ END